Inverse designed plasmonic metasurface with parts per billion optical hydrogen detection

Plasmonic sensors rely on optical resonances in metal nanoparticles and are typically limited by their broad spectral features. This constraint is particularly taxing for optical hydrogen sensors, in which hydrogen is absorbed inside optically-lossy Pd nanostructures and for which state-of-the-art detection limits are only at the low parts-per-million (ppm) range. Here, we overcome this limitation by inversely designing a plasmonic metasurface based on a periodic array of Pd nanoparticles. Guided by a particle swarm optimization algorithm, we numerically identify and experimentally demonstrate a sensor with an optimal balance between a narrow spectral linewidth and a large field enhancement inside the nanoparticles, enabling a measured hydrogen detection limit of 250 parts-per-billion (ppb). Our work significantly improves current plasmonic hydrogen sensor capabilities and, in a broader context, highlights the power of inverse design of plasmonic metasurfaces for ultrasensitive optical (gas) detection.


Ultrasensitive Hydrogen Sensors State-of-the-Art
In Table S1 we summarize the available literature on hydrogen sensors with different transduction principles (categorized following the definition given in Ref. 1 ) that report sensitivity down to the ppb level at room temperature. From the Table, excluding our results here, it is clear that until now, ppb sensitivity is only achieved by electrical (and one thermal conductivity) sensors.

Calculation of Rayleigh Anomalies in Periodic Nanoparticle Array
When light is incident on a 2D square array of subwavelength nanoparticle with pitch distance, a, photons can gain additional momentum in integer multiples of Rayleigh anomaly (RA) is associated with light diffracted parallel to the lattice surface, and it occurs when The wavelength of the RA occurs at the onset of the (i, j) diffraction order, above which free-space light diffraction is forbidden in the order. Here, is the incident angle.
Therefore, when the incident light has polarization of x direction, the wavelengths of (1,0) and (-1,0) orders of RAs are calculated by (±>,9) where is the refractive index of material at diffracted side and assuming the light is coming from air. Similarly, the wavelength of (0,±1) can be expressed as For simplicity, the wavelength or wavevector of higher orders of RAs also can be calculated from Eq. S3.

Angle Dispersion Extinction and Field Distribution of the Optimized Sensor Array
Supplementary Fig. 7. (a) Simulated wavelength-resolved angle dispersion extinction spectra of the optimized sensor array, showing the different RA orders (dashed lines) and the LSPR position of the corresponding single-particle counterpart (Supplementary Fig. 8). (b) Field distribution surrounding the nanoparticle (left) and at the close vicinity and inside of the nanoparticle (right) at three different excitation wavelengths corresponding to the extinction peaks (λ1-λ3). From the maps it is clear that relative field amplitude inside the nanoparticle excited at λ1 and λ2 are lower than the one at λ3. This again corroborates the nature of the peak, in which λ3 is dominated by the LSPR and thus is sensitive to the change from Pd to Pd hydride. Dashed lines outline the interfaces between glass/nanodisks/PMMA/air.

Single Particle Counterpart of the Optimized Sensor Array
Supplementary Fig. 8. Calculated Pd and PdH0.12 extinction spectra of a single Pd nanodisk with same geometrical parameters to the optimized array sensor obtained by PSO (d = 124 nm, h = 20 nm, a = 376 nm, tPMMA = 300 nm). Comparable peak shift as for the array sensor is found upon hydrogenation but with very broad FWHM. This combination results in a FoM of 0.07.

Extended Particle Swarm Optimization Method
If we assume the position of the Pd/PdH0.12 in the parameter space as a function of four variables, C ( ) = ( ( ), ℎ( ), ( ), HIIJ ( )), where d(n), h(n), a(n) and tPMMA(n) are the particle diameter, height and the pitch distance of the array, tPMMA is the thickness of PMMA, and n is the generation counter. During the optimizing process, the particle is subjected to three forces as it moves through the parameter spaces: (i) a frictional force that is proportional to the velocity, C ( ), where is the inertial weight; (ii) a spring force towards the individual best value of this particle, > > [ P C ( ) − C ( )], where > is the cognitive factor and > is a random number between 0 and 1; and (iii) a spring force towards the global best value of all the particles, ; ; [ ( ) − C ( )], where ; is the social factor and ; is a random number between 0 and 1. The velocity in the next generation can be obtained from the sum of these three forces, i.e., C ( + 1) = C ( ) + > > [ P C ( ) − C ( )] + ; ; [ ( ) − C ( )]. In the Lumerical solver, we used the default values of > = ; = 1.49 and linearly spaced values of between minimum 0.4 and maximum 0.9 for PSO simulations that have been verified to converge well in many test optimization for photonic design problems. The position of the particle in the next generation is then given by C ( + 1) = C ( ) + C ( + 1). Supplementary Fig. 9. PSO algorithm for one nanoparticle array in the PSO terminology updating during one generation. Supplementary Fig. 10. Simulated extinction spectra of Pd and PdH0.12 for Population 2 after 18 th generation of optimization. Clearly the PSO reaches a configuration where the two SLR peaks are close to each other and thus indistinguishable. Due to this condition, when the spectra change from Pd to Pd hydride, the program assigns a different peak, which then results in seemingly large peak shift (and thus falsely large FoM). To avoid this problem, we stopped our PSO simulation at the 15 th generation.  Supplementary Fig. 12. (a) Lorentzian function fitting (red dashed line) to the experimental optical spectra to extract λpeak. In our analysis, the fit is only applied within ±60 nm from the peak maximum (grey shaded area, following the method established in Ref. 17 ), where the peak is symmetric. (b) Zoomed-in version of (a) within the fit range. Clearly, the Lorentzian represents well the data and thus enables a good fit with R 2 > 0.99. (c) Lorentzian-fitted Δλpeak response of the best sensor (cf. Fig. 4a) in the first 30 min of operation used to derive the peak-to-peak readout noise, σ, of 0.01 nm. The dashed lines and gray-shaded areas denote the mean of the signal and ±σ from the mean, respectively. Supplementary Fig. 13. (a) Δλpeak response to stepwise random H2 concentration (250 to 0.25 ppm) in Ar carrier gas at room temperature. Inset: zoomed-in version of the sensor response to 250 ppb H2. (b) Measured Δλpeak as a function of H2 concentration derived from (a). The transparent symbols and gray dashed line are reproduced from Fig. 4b. The sensor's responses to these random H2 exposure are consistent with the descending one, and thus exemplifying the reproducibility of the sensor. Supplementary Fig. 14. (a) Δλpeak response to three consecutive cycles of 250 ppb H2 (grey areas). A reversible and reproducible sensor response to such low concentration of H2 is observed. (b) Average sensor signal to the three cycles of 250 ppb H2 exposure. An uncertainty of ~0.01 nm is recorded, which is in the same order of the sensor's signal noise.

Other Performance Aspects: Reproducibility, Speed, and Poisoning Resistance
To deduce the response and recovery times of the sensors we use the commonly used t90 and t10, respectively (Supplementary Fig. 15a-b). As plotted in Supplementary Fig. 15c-d, the response and recovery times of the sensor increase with the lowering H2 concentration. Such observation is inherent to Pd nanostructures as previously shown. 15,17,20 As a result, at the lowest H2 concentration of 250 ppb, the sensor's response time is in the order of 40 min. Interestingly, the response and recovery times of the control random array sensor are practically similar to the optimized one (Supplementary Fig. 15c-d). Such finding reveals that our method of increasing the sensor's sensitivity via periodic arrangement does not affect its sensing speed as it is mainly defined by the materials design. Thus, our sensitivity improvement method can be combined with other methods aimed to directly enhance the sensor's speed, for example by employing nanoparticles with reported faster kinetics than Pd (e.g. PdAu, 17,21 PdCo 15 and PdTa 22 alloys, with speed twice as fast compared to pure Pd) or by utilizing polymer coatings with higher kinetics-enhancements such as PTFE (twice as high as PMMA). 17 Supplementary Fig. 15. The definition of (a) response time as t90 and (b) recovery times as t10, which correspond to the time it takes to reach 90% and 10% of the normalized signal (with respect to signal during the exposure and in the absence of H2), respectively. (c) Response times and (d) recovery times of the optimized periodic array sensor and control random array sensor as function of H2 concentration. Data is extracted from Fig. 4a and c, respectively. The recovery and response times of both sensors are comparable and can practically be described with a single trend (the dashed lines), as established in ref. 17 . Supplementary Fig. 16. (a) Time-resolved Δλpeak response of sensor 1 pulse of 1000 ppm H2 followed by 5 pulses of 1000 ppm H2 + 500 ppm CO, and 1000 ppm H2 + 50 ppm NO2 in Ar. (b) Normalized sensor signal to the one obtained in 1000 ppm H2. The error bars denote the standard deviation from 5 cycles. The shaded area indicates the ±20% deviation limit from the normalized Δλpeak in 1000 ppm H2.  Fig. 5b) for Pd (light gray) and PdH0.12 (dark grey) nanodisk arrays. The spectra are basically identical to the ones of the sensor coated with 300 nm PMMA (cf. Fig. 3d).

Quasi-Random Array Control Sample: Structural, Optical and Noise Characterization
Supplementary Fig. 18. (a) SEM image of the control sensor array with quasi-random particle distribution. (b) Diameter distribution of the particles forming the sensor, showing an average of 122 nm, which is slightly smaller than the targeted diameter of 124 nm. (c) Radial distribution function (RDF) of the control sensor. The primary peak in the RDF (i.e., ~2.5⌀) indicates the average center-to-center distance between neighboring nanostructures. (d) Experimental extinction spectra of the quasi-random array.
Supplementary Fig. 19. (a) Lorentzian function fitting (red dashed line) within ±60 nm from the peak maximum (grey shaded area) to extract λpeak. (b) Zoomed-in version of (a) within the fit range. The Lorentzian represents the data well in the peak-maximum region and thus enables a good fit with R 2 > 0.97. (c) Lorentzian-fitted Δλpeak response of the quasi-random array control sensor (cf. Fig. 4c). The derived peak-to-peak noise, σcontrol, is 0.08 nm, much higher than that of the optimized regular array sensor. The dashed lines and gray-shaded areas denote the mean of the signal and ±σ from the mean, respectively.